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Random fields in nature often have, to a good approximation, Gaussian characteristics. We 
present the mathematical framework for a new and simple method for investigating the non-Gaussian 
contributions, based on counting the maxima and minima of a scalar field. We consider a random 
surface, whose height is given by a nonlinear function of a Gaussian field. We find that, as a result of 
the non-Gaussianity, the density of maxima and minima no longer match and calculate the relative 
imbalance between the two. Our approach allows to detect and quantify non-Gaussianities present 
in any random field that can be represented as the height of a smooth two-dimensional surface. 



A wide range of phenomena feature observables that 
can be regarded as random fields. The cosmic back- 
ground radiation [l| is a famous example, but the height 
profile of a growing surface , medical images of brain 
activity 0] and optical speckle patterns [1, Q also demon- 
strate this. 

In many cases, the fields can be approximated as Gaus- 
sian fields, meaning that they have certain properties 
which arc related to the Gaussian (or normal) distribu- 
tion. This is for example the case when the observable 
signal is averaged over a large scale, producing approxi- 
mately Gaussian statistics on account of the central limit 
theorem. The stochastic properties of such fields have 
already been the subject of several studies [H-Iioj: the 
density of maxima and minima for instance reflects the 
amount of field fiuctuations at short distances. 

Analytical investigations are often restricted to such 
Gaussian fields. However, phenomena described by non- 
linear laws produce non-Gaussian signals. Since these 
nonlinear effects are usually quite small, the resulting 
departures from Gaussianity can be tiny. Nevertheless, 
these non-Gaussianities can offer a key to understanding 
the interesting nonlinear processes behind the phenom- 
ena in question. 

If the non-Gaussianity is generated by microscopic 
nonlinear processes, then some indicator that is sensitive 
to short distances would be necessary to observe it. Mi- 
croscopic dynamics do not involve mixing between differ- 
ent regions, so the originally Gaussian field H{r) simply 
transforms in a local way, H{r) — ?> F{H{r)). Provided 
that this transformation is nonlinear, the new function 
will have non-Gaussian statistics. 

The standard approach to describing the statistics of 
a random field is to measure its correlation functions. In 
the case of a random scalar field h{x, y) with Gaussian 
statistics, its statistical properties are entirely encoded 
in its two-point correlation function {h{x,y)h{x' ,y')} (as 
a function of the distance between {x,y) and {x',y')). 
The higher-order correlation functions can be factorized 
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into two-point correlation functions, by Wick's theorem. 
A breakdown in these relationships is evidence that the 
field is not Gaussian. 

For example, to see that the field F{H) (as given 
above) is non-Gaussian, let us calculate its third-order 
correlation. Such a correlation would vanish (with re- 
spect to the mean) for a Gaussian variable. The third- 
order correlation function at equal points in space is the 
skewness {{F [H) - {F {H))f) . liF{H) = H + sH^, then 
the skewness is easily found to be 12£(i/^)^. Hence this 
field is non-Gaussian. 

Another way to measure the non-Gaussianity is to di- 
rectly measure the skewness. The equal-point correlation 
function should be most sensitive to non-Gaussianity be- 
cause F is local. This measurement also determines the 
value of e, so it gives some information about the dynam- 
ics of the nonlinear evolution. 

In this paper, we take a geometric approach to tackle 
this problem. We interpret the scalar field as the height 
of a surface (see fig. [J) and infer the statistical properties 
of the signa l by studying the stochastic topography of 
this surface [111 ■ Such an approach has already been the 
subject of both theoretical [a. ItI. [l^ - [l^ | and experimental 
studies Here, we focus on the statistical imbalance 
between peaks and troughs. A test of Gaussianity based 
on similar ideas has already been applied to the temper- 
ature fluctuations in the cosmic microwave background 

MM- 

We will focus on the difference between the densities 
of maxima and minima. This should also be sensitive 
to local statistics of the field, but it will be a measure- 
ment of the non-Gaussian properties in particular, since 
a Gaussian variable is always symmetric around its mean 
value. Wc will study signals of the form Fml{H) where 
the underlying field H is Gaussian and Fj^l is any nonlin- 
ear function, and we will find that the imbalance can be 
nonzero, illustrating this approach. Moreover, we show 
how large the imbalance is exactly in relation to the non- 
linear perturbation, which allows one to attack the re- 
verse problem: by measuring the difference in density 
between maxima and minima for a given near- Gaussian 
field, one can quantify the size of the non-Gaussian com- 
ponent. 
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Figure 1: A realization of a Gaussian field with periodic 
boundary conditions. 

The outline of this paper is as follows. In section U we 
review the properties of Gaussian fields and introduce 
the basic notions and notations that we will use. We 
then demonstrate how the imbalance between maxima 
and minima can be calculated in section HIl Section Hill is 
devoted to determining the key ingredient, namely the 
probability distribution for the values of minima in a 
Gaussian field. In section IIVI we arrive at the final re- 
sult, compare it with results from computer generated 
fields and point out the main features. Finally, section IVl 
provides an overview of our findings and their implica- 
tions. 



I. GAUSSIAN FIELDS 

The Gaussian distribution is the archetype of a contin- 
uous probability density. It is given by 

where fj, and a are the expectation value and standard 
deviation of the stochastic variable respectively. One of 
its special properties is that the sum of two indepen- 
dent stochastic variables, that adhere to this distribu- 
tion, is itself also a Gaussian variable, albeit of course 
with ^ = fii + ^2 and = crj 4- cr|. This property can 
be considered to be one of the components of the proof 
of the central limit theorem, which states that - under 
some very general conditions - the sum (or average) of 
a large number of independent stochastic variables ac- 
quires a Gaussian distribution, in the limit that the num- 
ber goes to infinity Because of this, many random 
processes can be well approximated using a Gaussian dis- 
tribution, e.g. the number of times a (fair) coin comes up 
heads when it is flipped a (large) number of times, or the 
amount of rain that falls at a certain spot during a year. 



A Gaussian random field is an extension of this princi- 
ple to two dimensions. For instance, one might consider 
the amount of rain that falls at different places through- 
out an area rather than a single spot. Upon adding to- 
gether all the contributions of all rain clouds during the 
course of a year, one obtains a random field. 

Formally, a field is a stochastic function H{r). The 
minimum requirement for a Gaussian field is that the 
probability distribution of H{'fo) at any point fo has to 
be described by a Gaussian. More generally, if we con- 
sider the values that the field attains at any number of 
points, ^i=ff(rl), ^=-^(^2), ^n=H{fn), the joint 
probability distribution has to be of the form 

p{^i, . . . ,^„) (X exp (^-^^A.j^.^j'^, (2) 

where Aij are constants. These constants give informa- 
tion about the relative values at different points (which 
would be useful for example if we wanted to know the 
distribution of the derivative of the field) . 

Any well-behaved Gaussian field can be decomposed 
into Fourier modes, resulting in the sum of an infinite 
number of wave functions 

i^(r) = V'o + ^^(fc)cos(fc-f-f 0,-). (3) 

k 

This shows how much of the fluctuations occur at each 
wavelength - for example, a surface of water might fluc- 
tuate with some random waves. If that is due to some 
external sound at a certain frequency, the Fourier trans- 
form will be strongest at the corresponding wavelength. 

This procedure may also be turned around - a Gaus- 
sian field may be generated by summing up a large num- 
ber of Fourier modes. We will now discuss a field that 
is generated in this way and try to understand how the 
statistics of the phase factors <j>j: reflect properties of the 
field, such as Gaussianity and translational invariance. 

The defining characteristic of a Gaussian field is now 
that the phases are random and completely uncor- 
related to each other. Already, by translational invari- 
ance, second order correlations between (jjk and (jik' are 
ruled out. If the phases arc completely independent, then 
at each individual point r, ipir) is the sum of an infi- 
nite number of independent random numbers between 
— 1 and 1 (as a result of the cosine), each weighted with 
a factor A{k). Thus, from the central limit theorem ip{f) 
is a Gaussian random variable. In contrast, in a non- 
Gaussian field the phases are correlated, i.e. the phases 
of different modes depend on each other. This mechanism 
is often called mode coupling. 

So far no statements have been made about the func- 
tion A{k): it has no infiuence on the Gaussianity (nor 
on the homogeneity) of ip. Indeed, this function is a 
free parameter, called the amplitude spectrum. While all 
Gaussian fields share some general properties, other more 
specific properties (such as the density of critical points, 
as we shall see) depend on this amplitude spectrum. For 
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example, when A(k) is large for vectors k with a small 
norm, the field ?/; is dominated by these waves with small 
wave vectors and hence large wave lengths, resulting in 
a more slowly varying ip as compared to a Gaussian field 
that is dominated by large wave vectors. 

There is one more condition that we will pose: next 
to being homogeneous, we will also only consider fields 
that arc isotropic, i.e. have rotational symmetry. This is 
achieved by requiring that A(k) depends on the magni- 
tude of k only, i.e. A{k) = A{k). 

In order to make a clear distinction between Gaus- 
sian and non-Gaussian, we will use H to indicate an 
(isotropic) Gaussian field and ip for any (homogeneous 
and isotropic) field. Later we will also use h, to indicate 
a perturbed Gaussian field. 

When we have a Gaussian variable x with a certain fi 
and cr, we can make a transformation to y = , which 
is then a standard Gaussian variable, having fi = and 
a — 1. This translation and rescaling has no effect on 
the overall properties of x and is introduced for conve- 
nience. We will apply a similar transformation, by set- 
ting (H) = and (i?^) = 1. The expectation values 
are obtained by integrating over all possible values of all 
random variables, which in this case, are the uniformly 
distributed phases: 

(,^,).(n/$)^^. (^) 

k 

For our earlier definition cq. ([3|) , the normalization trans- 
lates to Hq = and ^A{k)'^ = 1- This normalization 
is for the purpose of simplicity only and has no impact 
on our analysis. 

More details on these calculations, as well as additional 
properties of Gaussian fields and definitions, can be found 
in appendix 13 There we also demonstrate how the two- 
point correlation function can be derived from cq. ([3|). 
We also show how the higher-order correlation functions 
are related to the two-point ones. 

Testing whether these relations hold for a given field 
-0 can reveal whether is Gaussian or not. A more de- 
tailed analysis of the correlation functions can provide 
clues about the nature of the non-Gaussianity. 

Although correlation functions provide an excellent ap- 
proach from a purely mathematical point of view, deter- 
mining correlation functions for a given realization of a 
near-Gaussian field h may not always be practical, as it 
requires precise measurements of h in order to determine 
the correlation functions with a large enough precision. 

In this paper we consider a geometrical test for Gaus- 
sianity, which involves counting the number of maxima 
and minima. 



II. MAXIMA VERSUS MINIMA 

Due to symmetry, a Gaussian field H has as many min- 
ima as it has maxima. For a perturbed Gaussian field. 



like h = H + eH'^, this may no longer be the case. There- 
fore the difference in densities of maxima and minima can 
serve as an indication of non-Gaussianity. We shall now 
derive what this difference is in the generic case of a field 
given by /i(r) = FML{H{r)), where H \s a, Gaussian field 
and FjsiL is any (nonlinear) function (e.g. the identity 
plus a perturbation), which depends only on H(r), i.e. 
the original (unperturbed) value of the field at that same 
point. This scheme we will refer to as a local perturbation. 

Transforming the function with Jjvi does not move 
maxima and minima around, but it can interchange 
them, depending on the sign of F^j^ ~ dFj^L/dH at 
the point in question. To see this, note that maxima and 
minima, together with saddle points, are critical points. 
The critical points of h are given by 

= V/i(f) = ^Vi?(f) = F^^{H)WH{r). (5) 

We see that the critical points of H and h are the same 
points; however, the prefactor Fjyj^{H) may influence 
the type of critical point. The three types can be dis- 
tinguished by considering the second derivatives: saddle 
points have h^xhyy — h^y < 0, whereas for maxima and 
minima (together called extrema) this is positive. For 
maxima, unlike minima, we have h^x < - or hyy < 0. 

Consider a critical point ro and let z = H{r\)). The 
second derivatives of h at rp simply have an extra factor 
F'j^^{z) as compared to the second derivatives of H. This 
has no influence on the sign of hxxhyy — li^y, therefore, 
the saddle points (extrema) of H are also saddle points 
(extrema) of h. However, a maximum (minimum) of H 
is a minimum (maximum) of h when F'^j^(z) < 0. In 
order to determine how many extrema will undergo such 
a transformation, we need to know how often F^j^{z) < 
at such a point. 

Let g{z) be the probability density that a certain min- 
imum ro of H has the value H(r\j) = z. The probability 
P that a minimum of H becomes a maximum of h is then 

P = [ dzgiz). (6) 

For example, if we consider a square perturbation h = 
H + eH'^, for which F^]^{z) = 1 + 2ez, we have 

i_ 

P= / "dzg{z). (7) 

^ — oo 

Because of the symmetry of H, the maxima are dis- 
tributed according to g{—z). With that, we can simi- 
larly define a probability Q that a maximum becomes a 
minimum going from H to h. 

Let no be the density of minima (or maxima) of H. 
The density of minima (maxima) of H which are max- 
ima (minima) of h is then Png (Quq). We can quantify 
the resulting imbalance in maxima and minima in the 
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dimcnsionlcss parameter 

"-max - "rain (1 + P - Q)no - (1 - P + Q)nQ 



An : 



^niax I ^min 

P-Q 



2no 

<iz {g{z) - g{~z)). 



(8) 



Thus, if we can determine g{z), we can calculate the exact 
imbalance between the maxima and minima of h. 



III. DISTRIBUTION OF MINIMUM VALUES 

A. One dimension 

Let us first consider the probability distribution for 
minimum values of a Gaussian function on a line. We will 
then generalize to two dimensions, and afterward, discuss 
how these distributions depend on the power spectrum. 
We start with 



H{x) = ^ A{k) cos(fcx + 



(9) 



The minima are given by Hx{xq) = and H^xixo) > 0. 
We would thus like to know the probability density that 
H(xq) = z, given that Hx{xq) = and Hxx{xq) > 0: 

g{z) = p{H{Xmin) = z) 

= -p{H{xa) = zAHxixo) = A Hxx{xo) > 0). 
n 

(10) 

Here n = p{Hx{xo) = A Hxx{xo) > 0) can be identified 
as the density of the minima. 

We need to determine the joint probability distribution 
p{H(xo), Hx{xo), Hxx{xo)) - since H is homogeneous, p 
does not depend on ccq. 

Let us take a closer look at the first derivative 

Bx{xq) = A{k){-k) sm{kxQ + (j)k) 

k 

= ^fcA(fc)cos(fc.To + 0fc + i7r). (11) 

k 

We see that the expression for Hx still describes a Gaus- 
sian: the phases are simply increased by ^tt (modulo 2ti) 
and the spectrum has picked up a factor of k. The bot- 
tom line is that Hx{xo) is a Gaussian variable, and it is 
easy to confirm that the same goes for Hxx{xo) (or any 
derivative) . 

We thus have three Gaussian variables. The joint prob- 
ability distribution of a set of (correlated) Gaussian ran- 
dom variables is given by (compare cq. ^) 
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(27r)"/2VdrtC 



(12) 



Moreover, the coefficients C can be determined measur- 
ing the statistics of the field: it is the matrix of correla- 
tions 



(13) 



Let us calculate {H{x)Hxx{x)) as an example. Again, 
homogeneity allows us to set xq ~ for convenience. We 
then find 

{Hixo)Hxx{xo)) - {H{0)Hxxm 

A{k) COS (t)k J2 Mk')i-k'^) cos 

' k k' 

= A{k)A{k'){-k'^){cos ^k cos (j>k' ) 
fcfc' 

= ^A(fc)A(fc')(-fc")^4fc' 

kk' 

= X!-5^Wfc' = -^2. (14) 

k 

Here we made use of the moment K2 defined in eq. (|A7p . 

An even and an odd derivative of H are always uncor- 
related, e.g. 



(i/(0)if,(0)) =^^(fc)A(fc')(-fc')(cos</)fe sin(f>k') 

kk' 

= Y,Mkfi-k){cosq}kSmq}k)Skk' =0. (15) 

kk' 

This is because an even derivative features cosines while 
an odd derivative has sines, and their product averages 
to zero, as above. 

The final result is that for H, Hx and Hxx the corre- 
lations are 



C 



1 -K2 ' 

^2 

-K2 Ki 



(16) 



The determinant of C is K2(K4 — K2) and its inverse is 
1 



'K2K4, A'l^ 



/\ 4 - A'l 



(17) 







This gives 
p{H,Hx,Hxx) 



X exp 



2K, 



i2nr/^y^K2{K4-Ki) 

K^H^ + 2K2HHXX + H^x 



2{Ki - K. 



(18) 



The plan is now to set H = z and Hx = and integrate 
p over Hxx- However, one important factor still needs to 
be added. The probability wc have calculated is actually 
a probability density (since the probability that H'{xo) = 
and H{xo) = z exactly is zero), and it is not defined 
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with respect to the variables we need. It is defined by 
fixing a point Xq and determining the probability that 
Hx vanishes within a certain tolerance at that point: 

P(g(.To) £ [z,z + dz]AH4xo) € [0,dH']) 
dz dH' 

Instead, we actually want the probability that there is an 
exact critical point within a certain distance of xq- 



P 



3xm G [x[),xq +dx\ : 
H{x,n) e [z,z + dz] A Hx{x,n) = 
da; dz 

Over the range dx, dH' varies by 



dH' ^ 



dx 



dx = \Hxx\Ax. 



(19) 



In order to get the desired probability density with re- 
spect to x, we need to multiply our current probability 
density with 

The probability distribution for the minima is thus 
given by (see eq. (fTOj)) 



- [ dH,,, p{H = z,Hx^ 0, H,,) \H, 
n Jo 



(20) 



The prefactor, featuring the density of minima n, can 
be regarded as a normalization constant and is found by 
integrating g{z) over the entire z-range. This is easily 
accomplished by taking the expression above and first 
integrate over z, and only then over H,,. The result is 



/oo 
dzg{z)^l 
-OO 



1 



n=—^Iu/K2. (21) 

ZTT 



The integrand in eq. (PU)) is also Gaussian, but it is only 
integrated over for positive Hxx, resulting in 



1 - A 



27r 



■ exp 



z exp I 




A 



(22) 



2(1 -A) 

Here erfc is the complementary error function 



erfc(a;) = 



dte~ 



(23) 



which converges to 1 as x goes to — oo. The two param- 
eters K2 and K4 have been merged into a single dimen- 
sionless parameter 



A; 



Kl 



(0 < A < 1). 



Note that we set Kq = {H'^) = 1 for convenience. In the 
generic case Kq 7^ 1, we have A = K2/{KqK4). A proof 
that A < 1 is derived explicitly in the next section. 



B. Two dimensions 

In two dimensions, the procedure to calculate the dis- 
tribution of the minima is similar. The minima are 
defined by the conditions H, — Hy = (defining 
critical points), H^xHyy ~ H^y > (separating ex- 
trema from saddle points) and Hxx,Hyy > (distin- 
guishing minima from maxima). We thus need to find 
p{H, Hx, Hy, Hxx, Hyy, Hxy). Tfils Is stlll a Gaussian 
joint distribution function. 

We start again by determining the correlations, for ex- 
ample (again setting r = for convenience) 

(HxxHyy) = ^ A(fc)A(fc')fc^fc;2(cos0jcos0^,) 

kk' 

= Y,A{k)A{k')klk'^\5~^g = ^ \A{kfklkl 



kk' 



= — / dkdeU(k)k'^cos^esin'^e 

nOO 

= I / dkU{k)k^ = IK^. (25) 

In the third line we replaced the sum by an integral and 
performed it using polar coordinates. 

Remember from the one-dimensional case that the cor- 
relation of an even and an odd derivative is always zero, 
because in the calculation we encounter a product of a 
cosine and a sine, which integrated over the (random) 
phase yields zero. Based on the calculation method 
demonstrated above, we can make a more general state- 
ment: when the combined number of x-derivatives (y- 
derivatives) is odd, the integral over 9 (as above) features 
a cosine (sine) with an odd exponent; the integral over 
then gives zero. If we apply this rule to our six variables, 
we see that Hx, Hy and Hxy all have no "compatible 
match" in this respect; therefore, they are uncorrelated 
to all other variables. This allows us to factorize the joint 
probability distribution 



P{H, Hx,Hy, HxX,Hyy, Hxy) 

= p{Hx) p{Hy) p{Hxy) p{H, Hxx,Hy 



(26) 



The probability densities of the individual variables are 
straightforward, 

piHx)^^^o.,{-l-H!), (27a) 



p(F,) = ^i=exp(--i-i/2 



(27b) 



V^2 

p(i/..) = ^cxp(-i-<). (27c) 



(24) For H, Hxx and Hyy, we determine the correlation matrix 



1 -ii^2 



c 



IK4 




(28) 
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The determinant of C is jK^^Ki — i^Tf) and its inverse 



IS 



K2K4 

^2 



(29) 



K2K4, 2Kl - A 4 3^4 - 2Kl 
After some rearranging, eq. p2|) gives 

1 



X exp 



{K^H + + K2Hyy) 



(30) 



2iv:4 



Hi 



As in the one-dimensional case, we now have a proba- 
bihty density with respect to Hx and Hy, which we need 
to convert to one with respect to x and y. For that we 
need to multiply p with the Jacobian determinant 



d{Hx, Hy 



d{x,y) 



\HxxHyy 



(31) 



The probability distribution for the minima is thus given 
by 



g{z) = -p{Hx^Q)p{Hy = Q) 

dHxx'iHyydHxy P{H = Z,HxX,Hyy) 
Xp{Hxy)\HxxHyy-Hly\ 

dHxxd-HyydHxyP{z, HxX, Hyy) 

xp{Hxy)\HxxHyy-Hly\. (32) 



1 



mTK2 



The integrals must be taken over the volume for which 
HxxHyy — H^y > and Hxx,Hyy > 0, which forms the 
domain of the minima. These constraints and the inte- 
gration can be simplified by making the following change 
of variables. 



rCOS^ ^{HxX - Hyy), 

rsm9 = Hxy, 



dHxxdHyydHxy 



2\'-'-xx 1 '-'-yy) 

2rdrdsdd. 



(33a) 
(33b) 
(33c) 
(34) 



In terms of these new variables, we have HxxH, 



vv 



tt2 



r and the constraints of the volume arc given by 



< r < s. We get 

r27T POO PS 



1 



miK2 Jo Jo Jo 
X exp 



drdsdS ■ 



4r(s2 - r^) 



Kiz"^ + AK2SZ + As^ 



2{Ki - Kl) 



1^4 



' (35) 
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Figure 2: Histograms of the values of 10^ minima obtained 
from simulations, together with the distribution given by 
eq. P7)) . for (a) a disk spectrum (A = |); (b) a Gaussian 
spectrum (A — i). 



The density of the minima n can again readily be ob- 
tained by integrating over z: 



dzg{z) = 1 



(36) 



Note that this result matches the one obtained in [6|. 

After evaluating the double integral (taking care to 
integrate over r first), we obtain 



7 



27r(3- 2A) 
"3 



exp 



2(3 - 2A) 



' erfc 



A 



2(1-A)(3-2A) 



A(l - z^) exp (-hz^) crfc 



1 



A 



2(1 -A) 



(37) 



iv/3A(l-A)zexp^ 2(1 -A) 



The two parameters K2 and K4 have been merged into 
one as before, 



A; 



(0 < A < 1). 



(38) 



Again, when we set Kq = {H"^) 7^ 1, we get A = 
Ki/{KoIU). 

Let us prove that A < 1. After some rearranging, we 
see that this is equivalent to K0K4 — K2 > 0. We find 



K0K4 - A'l = / / dfcdfc' n(A:)n( 



{k)Il{k'){k"' -k-'k'-'). (39) 



Note that we could just as well replace fc'^ with k"^ (be- 
cause everything else is symmetric in k and /c'), and hence 
also with ^(fc^ + fc'^). If we do the latter, we can rewrite 



i(fc4 + fc'4)-fc2^'2 ^ 



(40) 



We see that this is positive, together with n(fc) and n(fc'), 
hence the integrand is positive and the integral too, which 
concludes the proof. 

We have compared eq. ([37]) with distributions obtained 
from computer-generated Gaussian fields - details about 
these numerical simulations and how the minima were 
identified can be found in appendix [BJ As can be seen in 
fig. [21 the agreement between cq. (|37]) and the numeric 
results is excellent. 

Let us take a closer look at eq. (jSTj • The two limits of 
A give results with interesting physical interpretations: 



lim 0(2) 

A— >0 



1 



-^^^-^/A. 



^/37r 



limg(z) = (l-sgnz)y — (e" 



1 



0(A), (41) 
)e-^^'. (42) 



The case A = occurs when K4 is unbounded (e.g. when 
n(/c) scales as fc-^). We see that the distribution is then 
an elementary Gaussian. A rough intuitive explanation 
for this is as follows. The key feature of this limit is that 
the maxima and minima arise from very rapid oscilla- 
tions that are superimposed on top of a slowly-varying 
field. In fact, if K4 is extremely large, the waves with 
a short wavelength (large \k\) have an amplitude that 
is small, but not negligible. They therefore create large 




Figure 3: Histogram of the values of 10^ minima obtained 
from simulations, together with the distribution given by 
eq. (|37p . for a ring spectrum (A = 1). No minima with a 
positive value of H were found. 



fluctuations in the gradient of the field and hence a lot of 
extrema; a fact that can also be seen from eq. (pS)) . Mean- 
while, the height of the surface at any point (including 
the abundant minima) is dominated by the waves with 
a large amplitude, which have long wavelengths (small 
|fc|). The location of the minima and the height of the 
surface are thus independent. Therefore, the distribution 
of the value of _ff at a minimum is the same as for any 
other point: Gaussian. 

Now we consider A = 1. From our proof that A < 1, 
it is not hard to see that this can only occur when 
n(fc) = 5{k — fco) for some constant fcg. This is called 
a ring spectrum, since the only occurring wave vectors 
are the ones with |fc| = fco, which describes a circle in 
fc-space. Inspecting eq. (|^^ we see that, due to the fac- 
tor (1 — sgnz), all minima have a negative value of _ff , as 
the simulations also show (see fig. [3]). The explanation 
is that height fields with a ring spectrum necessarily sat- 
isfy V^i? = —k^H - therefore, if H is positive, the mean 
curvature Hxx + Hyy < 0, so the point cannot be a min- 
imum. In other words, such Gaussian fields are random 
solutions to Helmholtz's equation - they could represent 
the height field of a large membrane resonating at a cer- 
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Figure 4: The skewness (71) and kurtosis (72) of the distri- 
bution (eq. (HTJ) as a function of A (see eqs. ((45)) and (|46)l '). 



tain frequency but with some randomness preventing a 
particular mode among the many at that frequency from 
stabihzing. 

While eq. (|37p appears quite complex, some of its pa- 
rameters have more transparent forms. The expectation 
value /X and standard deviation a for example are 



37r 



(43) 
(44) 



When looking at fig. [21 it appears that the distribution 
is itself almost Gaussian. This can be captured in the 
skewness 71 and kurtosis 72, 



71 



^3 _ 4^/2(64 - (18V3 - ll)7r) 

(37rA-i - (32 - (6\/3 - 2)77))^^^ 
3.46 



(9.42A-1- 5.63)3/2' 
--3 

4(-1536 + 32(18V^ - ll)7r + 9(2^3 - 9)77^ 



(45) 



M4 „ 
72 = — - 3 



(37rA-i - (32- (6V3-2)7r))^ 
2.68 



(9.42A-1 -5.63)2 



(46) 



Here ^„ is the n-th moment about the mean: ^„ = ((^ — 
(^))"). The skewness is a measure of the symmetry of a 
distribution around the mean, while the kurtosis gives an 
indication of its "peakiness". For a Gaussian distribution, 
both the skewness and the kurtosis are zero. They can 
therefore be considered as a measure of the Gaussianity 
of a distribution; note however that a distribution is not 
necessarily Gaussian if both parameters are zero. The 




(a) 




0.05 0.1 0.15 0.2 



(b) 



Figure 5: An for h — H + eH^ as a function of e, where H 
has a disk spectrum (A = |). The data points stem from 
simulations, the solid curve is eq. (|47p . The two graphs are 
for different ranges of e. 



two parameters are shown in fig. |4l Naturally, they both 
go to zero for A 0. 



IV. MAXIMA AND MINIMA IMBALANCE 

Now that we have obtained 17(2), we can calculate the 
relative imbalance between the densities of maxima and 
minima oi h ~ Fml{H), in accordance with eq. ([8]): 



An = 



dz {g{z) - g{-z)). 



(47) 
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0.05 0.1 0.15 0.2 
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(b) 



Figure 6: An for h — H + eH^ as a function of e, where H 
has a Gaussian spectrum (A = i). The data points stem from 
simulations, the solid curve is eq. H47p . The two graphs are 
for different ranges of e. 

The most basic example of a perturbed Gaussian for 
which we may expect An ^ is /i = iJ + eH^ . In this 
case, the domain of integration is [— oo,— We have 
compared eq. (|Tf)) with results from computer-generated 
fields, for two different spectra: in fig. [5] a so-called disk 
spectrum was used: 

A{kf^e{k-k^) K2n = ^ A = ^ (48) 

n + 1 4 

Fig. [S] features results for a Gaussian spectrum: 

A{kf - cxp(-fcV2fco) = 2"n!fc2" A = ^ 

(49) 



In both cases, we see an excellent agreement between 
the results from the simulations and our theoretical for- 
mula. 

In both figures, we see that Ati increases dramatically 
starting e ^ 0.15. This can be explained intuitively as 
follows: the balance in densities of maxima and minima is 
disturbed by extrema located below H = Since H is 
a standard Gaussian, such low values (i.e. large negative 
values) of H are exponentially rare. It is only when — ^ is 
in the order of —1 that a significant An can be expected. 
To get a rough estimate for the number of these extrema, 
we can just look at the density of points with H = 
(ignoring the requirement that they be minima does not 
change the exponential dependence) . This is e"^/^*'^ \ A 
more careful approximation (see appendix[C]) gives An ^ 

This argument also applies to the generic case h ~ 
H + sfNL{H), where /^r^ designates a perturbation and 

is a parameter controlling the size of the perturbation. 
Now e/]vL [H] needs to be in the order of 1 for An to be 
significantly nonzero. Thus measuring the imbalance be- 
tween maxima and minima does not give a very sensitive 
test of the type of non-Gaussianity that we have consid- 
ered here, in the limit of small e. However, eq. (|T7)) is a 
nonpcrturbative result that also holds for large e. 



V. CONCLUSIONS 

For a random field given by h{r) ~ Fn l{H (r)) , where 
H is a Gaussian field and Fnl any (nonlinear) function, 
we find that the densities of maxima and minima of h 
may differ. We have shown what the imbalance is as a 
function of the transformation Fjml and the power spec- 
trum of H. Our result is exact, and does not rely on 
perturbation theory, a nice feature since Fnl does not 
have to be small for our result to apply. This is con- 
firmed by our simulations. On the other hand, the im- 
balance between maxima and minima is exponentially 
small when e is small. Directly measuring the skewness 
at a given point, for example is much more sensitive to e 
when e <C 1. 

The simple reason is that, when H is of order one, 
h is a monotonic function of H, and hence it has the 
same numbers of maxima and minima. Only very large 
fluctuations in H can lead to an imbalance. Other types 
of non-Gaussian fields are more likely to have appreciable 
imbalances between maxima and minima. For example, 
the nonlinear evolution of a field, such as the height of 
a surface on which particles are accumulating could give 
rise to an imbalance between maxima and minima. The 
diffusion of the particles for instance might preferentially 
smooth out maxima. 
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Appendix A: Properties of Gaussian fields 

Let us start by calculating the mean and standard de- 
viation of a Gaussian field H, the equivalents of fi and 
cr of a Gaussian variable. This involves expectation val- 
ues, which are obtained by integrating over all possible 
values of all random variables, which in this case, are the 
uniformly distributed phases 



(Al) 



The mean is then simply 

(Hir)) ^(^Ho+J2 ^(^) cos(fc • r + ^g) 

k 

= Hn + Y^ A{k){cos{k ■ r + 0^)) 

k 

= Ho + Y, Mk) cos(fc • r + cj)j:) 

k 

= Ho. (A2) 
For the variance (standard deviation squared) we find 

{{H-{H)f) 

^ A{k) cos{k ■ r + 

k 

= ^A{k)A{k'){cos{k ■ r + (j)j:)cos{k' ■r + (j)g)). 

kk' 

(A3) 

Since the phases are uncorrelated, for fc 7^ fc' we find 

(cos(fc • f + <j)^) cos{k' ■ f+ 

= (cos{k-r + (/)^r))(cos(fc' • f+ <i)g)) = 0. (A4) 

Hence the term in the double sum can only be nonzero 
for k = k'. As a result we get 

{{H {H) f) = J2 A{kY{cos\k ■ r + 0,,)) 



(A5) 



For simplicity, we will set {H) = and (i/^) — 1, 
which translates to i/g = and J2k ^^i^)'^ — 1- 

While the vectors k in cq. ^ form a discrete set, usu- 
ally they are sufficiently finely spaced so that we can treat 
the amplitude spectrum A{k) as a continuous function 
defined over the positive reals. If we take our normaliza- 
tion condition, and replace the sum with an integral, we 
get 

Ak\a{kf = j j kdkde\a{kf 



/ dkTTka{kY = / dfcn(fc). 
/o Jo 



Here a{k) indicates the continuous spectrum equivalent 
to the discrete amplitudes A(k,). The newly introduced 
function n(A;) = Trka{k)'^ is the power spectrum of H. 

Some properties of a Gaussian field depend on the am- 
plitude spectrum. In many cases this dependence can be 
expressed in terms of the moments of the spectrum 

/•OO 

K,, = V iA(fc)2fc" = / dfcn(fc)fc". (A7) 
- Jo 



The normalization condition can be translated as Kq — 1. 



1. Two-point correlation function 

Correlation functions are often used to probe the Gaus- 
sianity of a given random field. This is because for Gaus- 
sian fields, they obey certain relations, as reviewed below. 
We will first calculate the two-point correlation function. 

The two-point correlation function C{fi,f2) of a field 
Tp is defined as 



(A8) 



When tjj is homogeneous and isotropic, C depends only 
on the distance between rl and 



C{R)^mr)iPif+R)), 



(A9) 



(A6) 



where r is any position and R is any vector of length 
R. For a Gaussian field H we find (if we set r = for 
convenience, which we are free to do) 

CiR) - {H{0)H{R)) 

A{k) cos(<?i^)) ( A{k) cos(fc -R + cf) 

k k 

= Y A{k)A{k'){cos{(j)f:) cos{k' ■ R + <^p)). (AlO) 

kk' 

Since the phases <j)^ arc uncorrelated, the correlation is 
automatically zero when k ^ fc', hence 

C{R) = Mk?{cos{4>j:) cos(fc • R + 0^)) 
fe 

= Y ^ikY{\ cos(A? • R + 2(j3j:) + \ cos(fc • R)). 

k 

(All) 

Since (/)^ is uniformly distributed, the expectation value 
of the first cosine is zero, and we are left with 

C{R) ^Y^Mkfcos{k-R) = J dk\a{kf cos[k- R). 

(A12) 

Wc thus find that the two-point correlation function of 
a Gaussian field is the Fourier transform of its (two- 
dimensional) power spectrum. Therefore, in essence, the 
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correlation function is as much a complete description of 
a Gaussian field as the power spectrum is. Also, by de- 
termining the correlation function and taking the inverse 
Fourier transform, one obtains the spectrum. 

The moments, defined before in terms of the power 
spectrum, can be related to the derivatives of the corre- 
lation function at R = 0. Because of symmetry, wc must 
have C{R) = C{—R). Hence C{R) is an even function 
and all its odd derivatives at zero vanish. To obtain the 
even derivatives, we must first eliminate the vector R in 
the equation above. We are free to choose its direction, 
so let us take R = Rx. We then get 



_d_ 

dR 



C(2")(0) = (-1)" Ak\a{kfk 



dk\a{kY cos{kxR) 



= / dfcia(/c)2(-l)"fc2"cos(fc^i?). 



(A13) 



= (-1)"/ / kdkd0^a{kfk'^''cos^"e 



(-i)"i^2„7^ fde cos'^e 

2tt J 
(2n~l)!! 



(AM) 



We thus find a one-to-one relation between the moments 
and the derivatives of the correlation function. The 
derivative of the correlation function C*-^"' (0) is related 
to roughness in the field itself. In fact, C*^^"''(0) is equal 
to (-l)"([i?(")(a;)]^), the fluctuations of the n-th deriva- 
tive (apart from a sign). 



2. Higher order correlation functions 

In general, the n-point correlation function is defined 
as the expectation value (V'(fl)V'(^^) • • ■ '0(''n)), as a func- 
tion of fi through r^. For a Gaussian field, this corre- 
lation function can be expressed in terms of two-point 
correlation functions, analogous to Wick's theorem. As a 
result, a non-Gaussian field can be recognized by check- 
ing whether this relation holds. In practice, this test can 
be applied to a single Gaussian field if it is homogeneous. 
In that case, the correlation function depends only on 
separations of the points — f\ through — . From a 
given homogeneous field "0, one obtains (a good approx- 
imation of) this correlation function by averaging over 
all (or a lot of) configurations with fixed spacings but 
translated to different rl's. 

The simplest case of the relationship is 

{HiH2H3Hi) = {HiH2){HsHi) + {HiH^){H2Hi) 
+ {HiHi){H2H3), 

(A15) 

where we introduced the notation Hi = H{fj) for short- 
ness. In general, correlations between an even number 



of variables with n > 2 can be reduced to the two-point 
correlations, while correlations between an odd number 
of variables always vanish. These properties follow from 
the definition of the Gaussian field: the n variables H{ri) 
are described by a correlated Gaussian distribution, and 
hence their correlation functions can be calculated ex- 
plicitly from Gaussian integrals. 

We shall now show how this characteristic relation 
comes about for our Fourier superposition. When we cal- 
culate the four-point correlation in the same way as we 
did for the two-point correlation, wc bring the brackets 
inside the (quadruple) sum, which gives us the term 



(C0s(fci • fi + C0s(fc2 • f2 + (^k^) 

X cos(fc3 • -f cos(fc4 • ''I + ^jj), 



(A16) 



which is summed for all combinations of fci through fc4. 
The first thing to note, is that whenever e.g. ki is not 
equal to any of the other ki, the correlation is automati- 
cally zero; this is because cos(fci • rl -|- 4>^^ ) is then inde- 
pendent of all other factors, can therefore be separated, 
and gives zero. Hence the correlation can only be nonzero 
if each fcj is equal to (at least) one other kt . We can dis- 
tinguish the cases fci = fc2, ^3 = ^4 and fci = ^3,^2 = ^4 
and fci = fc4, fc2 = k^. Let us focus on the first case; the 
sum of all these correlations gives 

^ A{kif A{k2,f {cos{ki ■ ri + cosik^ ■ ^ (jyj:^) 



ki ,^3 



X cos(fc3 • ri + ) cos(fc3 • rl + )). 



(A17) 



This can be split into 

^ A(fci)^(cos(fci • fl + 0gjcos(fci • r^2 + 

1 

^ A(fc3)^(cOs(fc3 • fi + (/«j^)cOs(fc3 • fl + (f)f:J) 



ki 
X 



^ {H{rl)H{ri)){H{ri)H{rl)). 



(A18) 



Applying the same to the other cases and adding them 
together precisely gives cq. (jA15[) . 

One may note that the case fci fc2 — ^3 = ^4 has 
not been treated correctly. However, since ki can take on 
an infinite number of values, and this case only provides 
one degree of freedom instead of the two we had for the 
other cases, these correlations only have an infinitesimal 
contribution. 

From this example, it is not hard to see that in general, 
an n-point correlation function can be factorized, that is, 
written as the sum of products of two-point correlations, 
where the sum features all possible ways in which the n 
variables can be paired up. 
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Figure 7: Identifying a critical point. The four black squares 
are grid points at which and Hy are known. At the red 
(blue) dots, = {Hy = 0) under the assumption that 
(Hy) is linear between the two grid points. The two contour 
lines Hx ~ and Hy = then intersect inside the square, 
indicating the existence of a critical point. 



Appendix B: Computer simulations 

In order to verify our theoretical results, we use a large 
number of computer-generated realizations of the Gaus- 
sian field H (typically a few thousand) - each with the 
same spectrum A{k) but random phases (j)^- We then 
apply the desired transformation Fnl and extract the 
desired statistics. 

The fields are defined on a square with periodic bound- 
ary conditions, i.e. H{x, y) = H(x + L,y) ~ H{x, y + L), 
in order to reduce finite size effects. This is accomplished 
by only adding together waves (as in cq. ([3])) with wave 
vectors k of which the x and y-components are multiples 
of 27r/L. 

The summation in eq. is restricted to wave vectors 
with a magnitude below a certain threshold k^ax- In 
order to minimize the potential effects of this cutoff, we 
choose spectra for which A(k) decays very quickly or is 
zero for large fc, such as the disk spectrum 



A{kf (xe{k~ko), 

and the Gaussian spectrum 

A(fc)2 cx exp(-fcV2fco)- 



(Bl) 



(B2) 



Finally, L is chosen in relation to k^ax such that (1) 
the sum in eq. ^ features at least a few hundred waves 
(recall that L influences this number via the periodic 
boundary conditions) and that (2) L is at least a few 
times 27r / ^/K2, which is a measure of the typical wave- 
length of the spectrum. An example of a Gaussian field 
generated in this way is shown in fig. [TJ 

The resulting formula for H is then evaluated at the 
grid points only. The distance between neighboring grid 
points is taken to be much smaller (by a factor of 50 
roughly) than the typical wavelength. Along with H it- 
self, we also calculate its first and second derivatives at 
these grid points. 

We use a very efficient approach to identifying criti- 
cal points and their type (maximum, minimum or saddle 
point). Every square of four neighboring grid points is 



considered. If or Hy has the same sign at all four 
points, we infer that it is not zero anywhere inside the 
square, which leads to the conclusion that the square does 
not contain a critical point. Otherwise, there would nec- 
essarily be at least two pairs of neighboring grid points 
with a different sign of H^. For each pair, it is as- 
sumed that Hx changes linearly between the two points, 
which allows to pinpoint two points along the edges of 
the square where ~ 0. The contour line = 
is then assumed to be a straight line between these two 
points. The same recipe is applied to Hy, after which 
it is determined whether the two contour lines crossed. 
The intersection (if present) is then a critical point. This 
idea is illustrated in fig. [T] 

It is also possible for all four neighboring points to 
have opposite signs of H^ (or Hy). This results in four 
points along the border of the square with H^ = 0, but 
without any information about which two pairs should 
be connected by a contour line. In combination with the 
two points with Hy = it can however be established 
what the parity of the number of intersections (i.e. critical 
points) is. We then simply assume this number to be or 
1 . This case is sufficiently rare (provided the grid is small 
enough) to not have a noticeable effect on the results. 

Once established whether the square under considera- 
tion contains a critical point, the type is determined by 
averaging the values of H^x, Hyy and H^y at the four 
grid points and evaluating the signs of H^xHyy — Hxy 
and Hxx + Hyy . 

Although this method clearly docs not always correctly 
determine the existence of a critical point or its type, it is 
not biased toward one outcome. Therefore, the mistakes 
that are made will get averaged out when statistics are 
taken over a large number of critical points. With regard 
to getting good statistics, the speed of this method is a 
big advantage. 

This method - along with the proper values of kmax, 
L and the grid size - was thoroughly tested on Gaussian 
fields (for which the statistical outcomes are known from 
theory) to verify its validity, before applying it to the 
non-Gaussian fields under investigation. 



Appendix C: Asymptotes for a very small 
non- Gaussianity 

In the limit where z is very large and negative, g(z) 
can be evaluated asymptotically. One can use the exact 
expression for g{z), but returning to the original integral 
eq. (|35|) gives more insight (and makes the calculations 
shorter) . The exponential weight in the integral is peaked 
at s = ^{Hxx + Hyy) = Therefore s almost cer- 

tainly becomes large and positive when — z is large, since 
the width of the distribution remains fixed. This allows 
us to extend the range of integration to all s's and to 
all r > 0, since the additional parts of the range have 
a very small weight. Then the integral can be worked 

I 2 

out exactly, giving g{z) « J -Az^e~T" apart from small 
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corrections, when -z is large and positive. We next sub- asymptotically (/^ z^e~^dx w Ae — ~ when A ^ <x>) 
stitute back in eq. (|S]). Noting that P dominates over Q • /\ pi \ 1^ 
and using integration by parts to evaluate the integral gives ^ ^ y 2ir e ^ 



[1] S. Dodelson, Modem Cosmology (Academic Press, 2003). 
[2] M. Kardar, G. Parisi, and Y. C. Zhang, Phys. Rev. Lett. 

56, 889 (1986). 
[3] K. J. Worsley, S. Marrett, P. Neelin, A. C. Vandal, K. J. 

Friston, and A. C. Evans, Human Brain Mapping 4, 58 

(1996). 

[4] F. Flossmann, K. O'HoUeran, M. R. Dennis, and M. J. 

Padgett, Phys. Rev. Lett. 100, 203902 (2008). 
[5] A. Weinrib and B. I. Halperin, Phys. Rev. B 26, 1362 

(1982). 

[6] M. S. Longuet-Higgins, Phil. Trans. R. Soc. Lond. A 250, 
157 (1957). 

[7] M. V. Berry and J. H. Hannay, J. Phys. A: Math. Gen. 

10, 1809 (1977). 
[8] G. Fohin, J. Phys. A: Math. Gen. 36, 1729 (2003). 



[9] G. Foltin, J. Phys. A: Math. Gen. 36, 4561 (2003). 
[10] I. Freund and M. Wilkinson, J. Opt. Soc. Am. A 15, 2892 
(1998). 

[11] R. Kamien, Rev. Mod. Phys. 74, 953 (2002). 
[12] M. R. Dennis, J. Phys. A: Math. Gen. 36, 6611 (2003). 
[13] M. Longuet-Higgins, Phil. Trans. R. Soc. Lond. A 249, 
321 (1957). 

[14] M. R. Dennis, Optics Letters 33, 2572 (2008). 

[15] A. F. Heavens and R. K. Sheth, Mon. Not. R. Astron. 

Soc. 310, 1062 (1999). 
[16] S. Gupta and A. F. Heavens, AIP Conf. Proc. 555, 337 

(2001). 

[17] N. G. van Kampen, Stochastic Processes in Physics and 
Chemistry (North-Holland Publishing company, 1981). 



